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DECAY BOUNDS FOR FUNCTIONS OF MATRICES WITH 
BANDED OR KRONECKER STRUCTURE 

MICHELE BENZI* AND VALERIA SIMONCINit 


Abstract. We present decay bounds for a broad class of Hermitian matrix functions where 
the matrix argument is banded or a Kronecker sum of banded matrices. Besides being significantly 
tighter than previous estimates, the new bounds closely capture the actual (non-monotonic) decay 
behavior of the entries of functions of matrices with Kronecker sum structure. We also discuss 
extensions to more general sparse matrices. 
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1. Introduction. The decay behavior of the entries of functions of banded and 
sparse matrices has attracted considerable interest over the years. It has been known 
for some time that if A is a banded Hermitian matrix and / is a smooth function 
with no singularities in a neighborhood of the spectrum of A, then the entries in /(A) 
usually exhibit rapid decay in magnitude away from the main diagonal. The decay 
rates are typically exponential, with even faster decay in the case of entire functions. 

The interest for the decay behavior of matrix functions stems largely from its 
importance for a number of applications including numerical analysis m m [m [n 
mi mi US] , harmonic analysis [H US] 132 , quantum chemistry [S] HUl |3S1 UJ , signal 
processing |341I42) . quantum information theory [nilTlim], multivariate statistics [T], 
queueing models [9] , control of large-scale dynamical systems |28j , quantum dynamics 
|24) . random matrix theory UD]j and others. The first case to be analyzed in detail 
was that of /(A) = A“^, see [TSl [T71 US 133] • In these papers one can find exponential 
decay bounds for the entries of the inverse of banded matrices. A related, but quite 
distinct line of research concerned the study of inverse-closed matrix algebras, where 
the decay behavior in the entries of a (usually infinite) matrix A is “inherited” by 
the entries of A~^. Here we mention [32]) where it was observed that a similar decay 
behavior occurs for the entries of /(A) = as well as [2 131 Uil USl 133] , among 

others. 

The study of the decay behavior for general analytic functions of banded matrices, 
including the important case of the matrix exponential, was initiated in [am] and 
continued for possibly non-normal matrices and general sparsity patterns in [7] ; further 
contributions in these directions include gmamiEi]- Collectively, these papers have 
largely elucidated the question of when one can expect exponential decay in the entries 
of /(A), in terms of conditions that the function / and the matrix A must satisfy. 
Some of these papers also address the important problem of when the rate of decay 
is asymptotically independent of the dimension n of the problem, a condition that 
allows, at least in principle, for the approximation of /(A) with a computational cost 
scaling linearly in n (see, e.g., 013 UQ]). 
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A limitation of these papers is that they provide decay bounds for the entries 
of f{A) that are often pessimistic and may not capture the correct, non-monotonic 
decay behavior actually observed in many situations of practical interest. A first step 
to address this issue was taken in where new bounds for the inverses of matrices 
that are Kronecker sums of banded matrices (a kind of structure of considerable 
importance in the numerical solution of PDE problems) were obtained; see also [89] 
for an early such analysis for a special class of matrices, and m for functions of 
multiband matrices. 

In this paper we build on the work in m to investigate the decay behavior 
in (Hermitian) matrix functions where the matrix is a Kronecker sum of banded 
matrices. We also present new bounds for functions of banded (more generally, sparse) 
Hermitian matrices. For certain broad classes of analytic functions that frequently 
arise in applications (including as special cases the resolvent, the inverse square root, 
and the exponential) we obtain improved decay bounds that capture much more 
closely the actual decay behavior of the matrix entries than previously published 
bounds. A significant difference with previous work in this area is that our bounds 
are expressed in integral form, and in order to apply the bounds to specific matrix 
functions it may be necessary to evaluate these integrals numerically. 

The paper is organized as follows. In section [5] we provide basic definitions and 
material from linear algebra and analysis utilized in the rest of the paper. In section[3] 
we briefly recall earlier work on decay bounds for matrix functions. New decay results 
for functions of banded matrices are given in section 01 Generalizations to more 
general sparse matrices are briefly discussed In section [SJ Functions of matrices with 
Kronecker sum structure are treated In section [SI Conclusive remarks are given in 
section jT] 

2. Preliminaries. In this section we give some basic definitions and background 
information on the types of matrices and functions considered in the paper. 

2.1. Banded matrices and Kronecker sums. We begin by recalling two stan¬ 
dard definitions. 

Definition 2.1. We say that a matrix M G is fi-banded if its entries M^j 

satisfy Mij = 0 for \i - j\ > j3. 

Definition 2.2. Let Mi, M 2 G C"^". We say that a matrix A G C" is the 
Kronecker sum of Mi and M 2 if 

A = Mi®M2-.= Mi®I + I®M2, ( 2 . 1 ) 

where I denotes the n x n identity matrix. 

In this paper we will be especially concerned with the case Mi = M2 = M, where 
M is /3-banded and Hermitian positive definite (HPD). In this case A is also HPD. 

The definition of Kronecker sum can easily be extended to three or more matrices. 
For instance, we can define 


A = Ml ® M 2 ® M^ ■.= Ml ® I ® I + I ® M 2 ® I + I ® I ® M^. 

The Kronecker sum of two matrices is well-behaved under matrix exponentiation. 
Indeed, the following relation holds (see, e.g., [29| Theorem 10.9]): 


exp(Mi © M2) = exp(Mi) © exp(M2). 


(2.2) 
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Similarly, the following matrix trigonometric identities hold for the matrix sine 
and cosine [IHl Theorem 12.2]: 

sin(Mi © M 2 ) = sin(Mi) © cos(M 2 ) + cos(Mi) © sin(M 2 ) (2.3) 

and 

cos(Mi © M 2 ) = cos(Mi) © cos(M 2 ) — sin(Mi) © sin(M 2 ). (2.4) 

As we will see, identity (12.21) will be useful in extending decay results for functions 
of banded matrices to functions of matrices with Kronecker sum structure. 

2.2. Classes of functions defined by integral transforms. We will be con¬ 
cerned with analytic functions of matrices. It is well known that if / is a function 
analytic in a domain C C containing the spectrum of a matrix A G C"^", then 

/(^) = (2.5) 

where i = \/—1 is the imaginary unit and T is any simple closed curve surrounding 
the eigenvalues of A and entirely contained in 17, oriented counterclockwise. 

Our main results concern certain analytic functions that can be represented as 
integral transforms of measures, in particular, strictly completely monotonie functions 
(associated with the Laplace-Stieltjes transform) and Markov functions (associated 
with the Cauchy-Stieltjes transform). Here we briefly review some basic properties 
of these functions and the relationship between the two classes. We begin with the 
following definition (see |44 ) ). 

Definition 2.3. Let f be defined in the interval {a,b) where —00 < a < b < + 00 . 
Then, f is said to be completely monotonic in (a, b) if 

(_l)'=y(fe)( 2 ;) > 0 for all a<x<b and all fc = 0,l,2,... 

Moreover, f is said to be strictly completely monotonic in (a, b) if 

(_l)'=y(fe)( 2 ;) > 0 for all a<x<b and all fc = 0,1,2,... 

Here denotes the fcth derivative of /, with f^^'i = f. It is shown in [33] that if 
/ is completely monotonic in (a, b), it can be extended to an analytic function in the 
open disk \z — h\ < b — a when b is finite. When b = + 00 , / is analytic in 5J(z) > a. 
Therefore, for each y G (a, b) we have that / is analytic in the open disk \z — y\ < R{y), 
where R{y) denotes the radius of convergence of the power series expansion of / about 
the point z = y. Clearly, R{y) > y — a for y G (a, b). 

In [5] Bernstein proved that a function / is completely monotonic in (0, 00 ) if and 
only if / is the Laplace-Stieltjes transform of q;(t); 

pOO 

f{x) = / e"™da(T), (2.6) 

Jo 

where a{T) is nondecreasing and the integral in (12.61) converges for all x > 0. Moreover, 
under the same assumptions / can be extended to an analytic function on the positive 
half-plane 5R(z) > 0. A refinement of this result (see [20]) states that / is strictly 
completely monotonic in (0, 00 ) if it is completely monotonic there and moreover the 
function Q!(r) has at least one positive point of increase, that is, there exists a tq > 0 
such that a(ro + (5) > Q;(ro) for any i5 > 0. 

Prominent examples of strictly completely monotonic functions include (see |43j 1: 
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1 - fi{x) = 1 /x = /g°° e ^'^dai{T) for cc > 0 , where ai(r) = r for r > 0 . 

2 . f 2 {x) = e~^ = e~^'^da 2 {T) for x > 0 , where a 2 (T) = 0 for 0 < r < 1 and 

a 2 (T) = 1 for T > 1 . 

3. fs^x) = (1 — e~^)/x = e““'^dQ! 3 (T) for x > 0, where a 3 (T) = r for 

0 < T < 1, and a 3 (r) = 1 for r > 1 . 

Other examples include the functions x~'^ (for any cr > 0), log(l + 1/x) and 
exp(l/a:), all strictly completely monotonic on (0, oo). Also, products and positive 
linear combinations of strictly completely monotonic functions are strictly completely 
monotonic, as one can readily check. 

A closely related class of functions is given by the Cauchy-Stieltjes (or Markov- 
type) functions, which can be written as 

/W=/hM, ,,c\r, (2,7) 

Jr z-oj 

where 7 is a (complex) measure supported on a closed set F C C and the integral is 
absolutely convergent. In this paper we are especially interested in the special case 
r = (— 00 ,0] so that 


fix) 



X G C \ (— 00 ,0], 


where 7 is now a (possibly signed) real measure. The following functions, which occur 
in various applications (see, e.g., [27] and references therein), fall into this class: 


0 


1 


1 


zdw, 


, Z — W TTyf—OJ 

\ rO 


log(l z) 


' —00 
/*-! 


1 Shl{ty/—Uj) 
Z — U! —TTUJ 

^ ^ dw. 


dw. 


' — LU (—w) 


The two classes of functions just introduced overlap. Indeed, it is easy to see 
(e.g., [3B]) that functions of the form 


fix) = f 

^0 


dA^(5) 

X + Oj' 


with /X a positive measure, are strictly completely monotonic on ( 0 , 00 ); but every 
such function can also be written in the form 


fix) 



7(w) 


-fii-uj), 


and therefore it is a Cauchy-Stieltjes function. We note that f{x) = exp(—x) is an 
example of a function that is strictly completely monotonic but not a Cauchy-Stieltjes 
function. 

In the rest of the paper, the term Laplace-Stieltjes function will be used to denote 
a function that is strictly completely monotonic on ( 0 , 00 ). 
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3. Previous work. In this section we briefly review some previous decay results 
from the literature. 

Given a n x n Hermitian positive definite /3-banded matrix M, it was shown in 
[TT] that 

\{M-%\<Cq^ (3.1) 

for all i,j = 1 ,... ,n, where q = {y/n — l)/{y/K + 1), k is the spectral condition 
number of M, C = max{l/Ai„in(.^)) C}) and (7 = (1 + . In this 

bound the diagonal elements of M are assumed not to be greater than one, which 
can always be satisfied by dividing M by its largest diagonal entry, after which the 
bound m will have to be multiplied by its reciprocal. The bound is known to 
be sharp, in the sense that it is attained for a certain tridiagonal Toeplitz matrix. 
We mention that is also valid for infinite and bi-infinite matrices as long as 

they have finite condition number, i.e., both M and M~^ are bounded. Using the 
identity M~^ = {M*M)~^M*, simple decay bounds were also obtained in [T7] for 
non-Hermitian matrices. 

Similarly, if M is /3-banded and Hermitian and / is analytic on a region of the 
complex plane containing the spectrum cr(M) of M, then there exist positive constants 
C and q < I such that 

\{f{M)U<Cq^, (3.2) 

where C and q can be expressed in terms of the parameter of a certain ellipse sur¬ 
rounding a{M) and of the maximum modulus of / on this ellipse; see [^. The bound 
(1^ . in general, is not sharp; in fact, since there are infinitely many ellipses contain¬ 
ing a{M) in their interior and such that / is analytic inside the ellipse and continuous 
on it, one should think of (13.21) as a parametric family of bounds rather than a single 
bound. By tuning the parameter of the ellipse one can obtain different bounds, usu¬ 
ally involving a trade-off between the values of C and q. This result was extended in 
[7] to the case where M is a sparse matrix with a general sparsity pattern, using the 
graph distance instead of the distance from the main diagonal; see also [nisa and 
section [5] below. Similar bounds for analytic functions of non-Hermitian matrices can 
be found in [H |7] . 

Practically all of the above results consist of exponential decay bounds on the 
magnitude of the entries of f{M). However, for entire functions the actual decay is 
typically superexponential, rather than exponential. Such bounds have been obtained 
by Iserles for the exponential of a tridiagonal matrix in m- This paper also presents 
superexponential decay bounds for the exponential of banded matrices, but the bounds 
only apply at sufficiently large distances from the main diagonal. None of these bounds 
require M to be Hermitian. Superexponential decay bounds for the exponential of 
certain infinite tridiagonal skew-Hermitian matrices arising in quantum mechanical 
computations have been recently obtained in |41| . 

4. Decay estimates for functions of a banded matrix. In this section we 
present new decay bounds for functions of matrices f{M) where M is a banded, 
Hermitian and positive definite. First, we make use of an important result from |30) 
to obtain decay bounds for the entries of the exponential of a banded, Hermitian, 
positive semidefinite matrix M. This result will then be used to obtain bounds or 
estimates on the entries of /(M), where / is strictly completely monotonic. In a 
similar manner, we will obtain bounds or estimates on the entries of f{M) where / 
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is a Markov function by making use of the classical bounds of Demko et al. [T7] for 
the entries of the inverses of banded positive definite matrices. 

In section ini we will use these results to obtain bounds for matrix functions f{A), 
where M is a Kronecker sum of banded matrices and / belongs to one of the two 
above-mentioned classes of functions. 

4.1. The exponential of a banded Hermitian matrix. We first recall (with 
a slightly different notation) an important result due to Hochbruck and Lubich [50] . 
Here the m columns of Vm G form an orthonormal basis for the Krylov subspace 

Kjn{M, v) = span{?;, Mv, ..., with ||u|| = 1, and Hm = V^MVm- 

Theorem 4.1. Let M be a Hermitian positive semidefinite matrix with eigenval¬ 
ues in the interval [0,4p]. Then the error in the Arnoldi approximation of exp{TM)v 
with ||u|| = 1, namely Sm '■= || exp(—rM)u — exp(—||, is bounded in the 
following ways: 

i) Sm < 10 exp(—m^/(5pT)), for pr > 1 and y/Apr < m < 2pT 
ii) Em < 10(pr)“i exp(-pT) (^)™ for m > 2pr. 

With this result we can establish bounds for the entries of the exponential of a 
banded Hermitian matrix. 

Theorem 4.2. Let M be as in Theorem \ f.l\ Assume in addition that M is 
fd-handed. Then, with the notation of Theorem \4.1\ and for k ^ t: 
i) For pr > 1 and y/Apr < \k — t\/fi < 2pr, 


\{exp{-TM))kt\ < 10exp 


5pr J ’ 


ii) For \k — t\/jd >2pT, 


|(exp(-rM))fet| < 10 


exp(-pr) 

pr 



Proof. We first note that an element of the Krylov subspace is a 

polynomial in M times a vector, so that 1/^ exp(—= Pm-iiTM)v for some 
polynomial Pm-i of degree at most m — 1. Because M is Hermitian and /3-banded, 
the matrix pm-iirM) is at most (m — l)/3-banded. 

Let now k,t with k ^ t he fixed, and write \k — t\ = (m — l)/3 -|- s for some 
TO > 1 and s = in particular, we see that {pm-i(TM))kt = 0 , moreover 

\k — t\/fd < TO. Consider first case ii). If to > 2pr, for u = et we obtain 


|(exp(-rM))fet| = |(exp(-rM))fet - {pm-i{TM))kt\ 

= \sl{exp{-TM)et - Pm-i{'rM)et)\ 

< ||exp(-TM)et-pm-i(TAL)et|| 

< 10(prr> exp(-pr) 

where in the last inequality Theorem 14.11 was used for m\k — t\/fd > 2pr. An analogous 
result is obtained for to in the finite interval, so as to verify i). □ 

As remarked in [30], the restriction to positive semidefinite M leads to no loss 
of generality, since a shift from M to M 5L entails a change by a factor in the 
quantities of interest. 
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Fig. 4.1. Example \4.3\ Bounds for \ ex.p{—T{M — X^iYi^))\:,t, t = 127, using Theorem \4.2\ M of 
size n = 200 and r = A. Left: M = tridiag(—1, 4, —1). Right: M = pentadiag(—0.5, —1, 4, —1, —0.5). 
Logarithmic scale. 


We also notice that in addition to Theorem 14.11 other asymptotic bounds exist for 
estimating the error in the exponential function with Krylov subspace approximation; 
see, e.g., HUH]. An advantage of Theorem 14.11 is that it provides explicit upper 
bounds, which can then be easily used for our purposes. 

Example 4.3. Figure ITTI shows the behavior of the bound in Theorem 14.21 for 
two typical matrices. The plot on the left refers to the tridiagonal matrix M = 
tridiag(—1,4, —1) (/3 = 1) of size n = 200, with r = 4, so that rp Ri 3.9995. The tth 
column with t = 127 is reported, and only the values above 10“®° are shown. The plot 
on the right refers to the pentadiagonal matrix M = pentadiag(—0.5, —1,4, —1, —0.5) 
(/3 = 2) of size n = 200, with r = 4, so that rp r; 4.4989. The same column t = 127 is 
shown. Note the superexponential decay behavior. In both cases, the estimate seems 
to be rather sharp. 


4.2. Bounds for Laplace—Stieltjes functions. By exploiting the connection 
between the exponential function and Laplace-Stieltjes functions, we can apply The¬ 
orem 14.21 to obtain bounds or estimates for the entries of Laplace-Stieltjes matrix 
functions. 

Theorem 4.4. Let M = M*P-banded and positive definite, and let M = 
M — Ainind, with the spectrum of M contained in [0,4p]. Assume f is a Laplace- 
Stieltjes function, so that it can be written in the form f{x) = Jf‘e~^'^da{T). Then, 
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with the notation and assumptions of Theorem \4.S\ and for \k — t\/fi > 2: 

POC 

\f{M)\k,t< / exp(-AininT)|(exp(-rM))fc,i|da(r) 


JO 

< 10 /'‘^»p(-A„,.r) 




/o 


pr 


\k-t\ 

0 


da{T) (4-1) 


+ 10 


/ 4p/3H 

l \k-t\ 

~2W 


exp(-Aminr) exp - 


i\k-t\/py 


bpr 


da(T) 


+ / exp(-Ai„inT)(exp(-rM))fc tda(r) =/ + // + III. 

4^ 


In general, these integrals may have to be evaluated numerically. We observe that 
in the above bound, the last term (III) does not significantly contribute provided that 
\k — t\ is sufficiently large while p and /3 are not too large. 

As an illustration, consider the function f{x) = Ijy/x. For this function we have 
Q;(r) = -^r(—^ + 1) = x/tt/t with r G (0,oo). We have 


I + II = 10^/n 


. 1 ^_IL / •. 

2 p^ exp(—An 


it) exp(—pr) j epr 


lk-t\ 




pr 


\k-t\ 


dr 


\k-t\ 


+ 10v^ 


exp(-AminT) 


\k-t\ 

2pfi 


-y/r 


exp 


5pT 


dr. 


while 


III < 



exp(-Ai„inT) , 
- -j= -dr. 


Figure 1121 shows two typical bounds for the entries of M~^ for the same matrices 
M considered in Example 14.31 The integrals I and II and the one appearing in 
the upper bound for III have been evaluated accurately using the built-in Matlab 
function quad. Note that the decay is now exponential. In both cases, the decay is 
satisfactorily captured. 

As yet another example, consider the entire function f{x) = (1 — exp(—x))/a; 
for a; G [0,1], which is a Laplace-Stieltjes function with da(r) = dr (see section 221) • 
Starting from 63 ) we can determine new terms I, II, and estimate III as it was done 
for the inverse square root. Due to the small interval size, the first term I accounts 
for the whole bound for most choices of k, t. For the same two matrices used above, 
the actual (superexponential) decay and its approximation are reported in Figure 11731 

We remark that for the validity of Theorem 14.41 we cannot relax the assumption 
that M be positive definite. This makes sense since we are considering functions of M 
that are defined on (0,oo). If M is not positive definite but / happens to be defined 
on a larger interval containing the spectrum of M, for instance on all of R, it may 
still be possible, in some cases, to obtain bounds for f{M) from the corresponding 
bounds on f{M + 61) where the shifted matrix M + SI is positive definite. 

Remark 4.5. We observe that if f{M + if I) is well defined for <7 G M, then the 
estimate also holds for \f{M + ifl)\k.t, since \ exp(iC)| = 1. 
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Fig. 4.2. Estimates for \M t = 127, using I+II and the upper bound for III. Size 

n = 200, Log-scale. Left: M = tridiag(— 1, 4, —1). Right: M = pentadiag(—0.5, — 1, 4, — 1, —0.5). 



Fig. 4.3. Estimates for \M ^(J—exp(—t = 127, using I-h II and the upper bound for III. 
size n = 200, Log-scale. Left: M = tridiag(—1, 4, —1). Right: M = pentadiag(—0.5, —1, 4, —1, —0.5). 


4.3. Bounds for Cauchy—Stieltjes functions. Bounds for the entries of f{M), 
where / is a Cauchy-Stieltjes function and M = M* is positive definite, can be ob¬ 
tained in a similar manner, with the bound (EH) of Demko et al. HZ] replacing the 
bounds on exp(— tM) from Theorem 14.21 

For a given w e F = (oo,0), let k = k{uj) = (Amax - a;)/(Amin - w), g = q{u}) = 
{^/k - 1 )/(\/k -b 1), C = C{-uj) = max{l/(Ai„in - with Co = Co(-a;) = 

(1 -b j (2(Amax — u})). We immediately obtain the following result. 

Theorem 4.6. Let M = M* be positive definite and let f be a Cauchy-Stieltjes 
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function. Then for all k and t we have 


\f{M)kt\< 


fc —1| 


C(uj)q(uj)~^ dw 


(4.2) 


For specific functions we can be more explicit, and provide more insightful upper 
bounds by evaluating or bounding the integral on the right-hand side of (14.21) . As 
an example, let us consider again f{x) = x~^, which happens to be both a Laplace- 
Stieltjes and a Cauchy-Stieltjes function. In this case we find the bound 

IA47I < -(C(0) + c,) (' , ( 43 ) 

TT \ Y ^max H" V Amin / 

where C 2 = max |l, (1 -f |. Indeed, for the given function and upon sub¬ 

stituting T = —oj, (I4.2[l becomes 


1^.71 < 


C'(r) 


aA max \/Amin + T 

vA max T \/Amin “t“ 'T' 






dr 


< 


1 f ^ max yA n 

TT 

max -I- vA n 


\k-t\ 

13 


C{T)-^dT. 

vA 


(4.4) 

(4.5) 


Let (()(r) be the integrand function. We split the integral as (/)(T)dr = (f>(T)dT + 
(/)(T)dT. For the first integral, we observe that C'(t) < (7(0), so that 


f C { T )^ dT < C { 0 ) [ ^dT = 2(7(0). 
Jo 7 Jo 7 


For the second integral, we observe that C{t) < € 2 ^ where C 2 = max{l,(l -|- 
\/k(0))^A/2}, so that 



(7(r)^dT < C 2 

7 



^dr = 2(72. 

Ty/T 


Collecting all results the final upper bound (14.31) follows. 

We note that for this particular matrix function, using the approach just presented 
results in much more explicit bounds than those obtained earlier using the Laplace- 
Stieltjes representation, which required the numerical evaluation of three integrals. 
Also, since the bound dSH) is known to be sharp (see [H]), it is to be expected that 
the bounds (14.31) will be generally better than those obtained in the previous section. 
Figure 14.41 shows the accuracy of the bounds in (14.31) for the same matrices as in 
Figure 14.21 where the Laplace-Stieltjes bounds were used. For both matrices, the 
quality of the Cauchy-Stieltjes bound is clearly superior. 

We conclude this section with a discussion on decay bounds for functions of M — 
ifl, where C G R. These estimates may be useful when the integral is over a complex 
curve. We first recall a result of Freund for {M — To this end, we let again 

Amin,Aniax be the extreme eigenvalues of M (assumed to be HPD), and we let Ai = 
Amin fC? -^2 — Amax ^C* 

Proposition 4.7. Assume M is Hermitian positive definite and ft-banded. Let 
R > 1 be defined as R = a + y/a'^ — 1, with a = (|Ai| + |A 2 |)/|A 2 — Ai|. Then for 
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Fig. 4.4. Estimates for \M t = 127, using 114-3]) , size n = 200, Log-scale. Left: 

M = tridiag(—1,4, —1). Right: M = pentadiag(—0.5, —1, 4, —1, —0.5). 


k^t, 

\t-k\ 

|(M - < CO ( 1 ) ^ w.th C(0 = 

With this bound, we can modify (14.21) so as to handle more general matrices as 
follows. Once again, we let A„iin, Amax be the extreme eigenvalues of M, and now we 
let Ai = Amin — iC, — uj 1 \2 = Amax — K ~ Q! and R are defined accordingly. 

i/(M-</)ia< r c 

a —oo 

Since R = R{C, w) is defined in terms of spectral information of the shifted matrix 
M — i(I — ojI, we also obtain C = C((^,uj) = 

^ / |A„ax-Amin| -1) 

5. Extensions to more general sparse matrices. Although all our main 
results so far have been stated for matrices that are banded, it is possible to extend 
the previous bounds to functions of matrices with general sparsity patterns. 

Following the approach in [14] and |7], let G = {V,E) be the undirected graph 
describing the nonzero pattern of M. Here H is a set of n vertices (one for each 
row/column of M) and E is a set of edges. The set E GV xV \s defined as follows: 
there is an edge {i,j) G E if and only if My ^ 0 (equivalently, Mji ^ 0 since 
M = M*). Given any two nodes i and j in V, a path of length k between i and j 
is a sequence of nodes iq = = j such that G E for all 

£ = 0, l,...,/c — 1 and ii im for I ^ m. If G is connected (equivalently, if M is 
irreducible, which we will assume to be the case), then there exists a path between 
any two nodes i,j € V. The geodesic distance d(i,j) between two nodes i,j & G is 
then the length of the shortest path joining i and j. With this distance, (G, d) is a 
metric space. 

We can then extend every one of the bounds seen so far for banded M to a general 
sparse matrix M = M* simply by systematically replacing the quantity by the 


d7(w), k ^ t. 


(4.6) 
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geodesic distance d{k, t). Hence, the decay in the entries of f{M) is to be understood 
in terms of distance from the nonzero pattern of M, rather than away from the main 
diagonal. We refer again to [7] for details. We note that this extension easily carries 
over to the bounds presented in the following section. 

Finally, we observe that all the results in this paper apply to the case where M is 
an infinite matrix with bounded spectrum, provided that / has no singularities on an 
open neighborhood of the spectral interval [Amin, Amax]- This implies that our bounds 
apply to all the n x n principal submatrices (“finite sections”) of such matrices, and 
that the bounds are uniform in n as n —>■ oo. 

6. Estimates for functions of Kronecker sums of matrices. The decay 
pattern for matrices with Kronecker structure has a rich structure. In addition to 
a decay away from the diagonal, which depends on the matrix bandwidth, a “local” 
decay can be observed within the bandwidth; see Figure lOI This particular pattern 
was described for f{x) = x~^ in [TT]; here we largely expand on the class of functions 
for which the phenomenon can be described. 



Fig. 6.1. Three-dimensional decay plots for [f{A)]ij where A. is the 5-point finite difference 
discretization of the negative Laplacian on the unit square ora a 10 X 10 uniform grid with zero 
Dirichlet boundary conditions. Left: f{A-) = exp(—5yt). Right: f{A-) = ■ 


Some matrix functions enjoy properties that make their application to Kronecker 
sums of matrices particularly simple. This is the case for instance of the exponential 
and certain trigonometric functions like sin(a;) and cos(a:). For these, bounds for their 
entries can be directly obtained from the estimates of the previous sections. 

6.1. The exponential function. Recall the relation (12.2L which implies that 

exp(—T^) = exp(—rM) 0 exp(— tM), r G K (6.1) 

when A = M ® 1 +1 ® M. Here and in the following, a lexicographic ordering of the 
entries will be used, so that each row or column index k oi A corresponds to the pair 
k = (fci, ^ 2 ) in the two-dimensional Cartesian grid. Furthermore, for any fixed values 
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of t,p,I3 > 0, define 




10 exp 


(- 


5pr 


IQ exp(-pr) 

pr 




for < 2pT, 

for > 2pT. 


( 6 . 2 ) 


Note that j) is only defined for \i — j\ > y4pr/3. With these notations, the 
following bounds can be obtained. 

Theorem 6.1. Let M he Hermitian and positive semidefinite with bandwidth /3 
and spectrum contained in [0,4p], and let A = I ® M + M ® I. Then 


(exp(-Ty4))fct = (exp(-TM))fcjti(exp(-rM))fe2t2. 


Therefore, for t > 0 


|(exp(-T^))fct| < 4>(fcl,tl)$(fc2,t2) 


for all t = (^ 1 ,^ 2 ) and k = (fci, ^ 2 ) with min{|ti — fci|, |t 2 — fe|} > 

Proof Using (1^ we obtain 

exp(—r^)et = exp(— tM) 0 exp(—rM)et. 

Let Et^t 2 be the n x n matrix such that e* = vec{Et^t 2 ) G i in particular 
Etj_t 2 = with eti,et 2 G M". Then 


ef exp{—TM) (g) exp(—rM)e( = e^vec{exp{—TM)Et^t 2 exp(—rM)*) 

= e^vec(exp(—rM)e(je^ exp(—rM)*) 

exp(-rM)eti(e| exp(-rM)*)ei 
exp{-TM)et^{ef^ exp(-rM)*)e 2 


- p'^ 


exp(-rM)eti(ef^ exp(-rM)*)e„ 
= ef, exp(-rM)eti(e^^ exp(-rM)*)efeJ, 


which proves the first relation for M Hermitian. For the bound, it is sufficient to use 
(j6.2() to obtain the desired conclusion. □ 

The result can be easily generalized to a broader class of matrices. 

Corollary 6.2. Let A = L® Mi + M 2 0 1 with Mi and M 2 having bandwidth jdi 
and P 2 , respectively. Also, let the spectrum of Mi be contained in the interval [0,4pi] 
and that of M 2 in the interval [0,4p2]j with pi, P 2 > 0. Then for t = (^ 1 ,^ 2 ) and 
k = {ki,k 2 ), with \ti - k(\ > £ = 1,2, 


{exp{-TA))kt = (exp(-rMi))fcjti(exp(-rM2))t2fe2- 


Therefore, 


|(exp(-ry4))fet| < $i(fci,ti)$ 2 (fc 2 ,t 2 ) 

where ^k{i,j) defined as 4>(i, j) in hh.Al with pi, (3i replacing p, 13. 

Generalization to the case of Kronecker sums of more than two matrices is rel¬ 
atively straightforward. Consider for example the case of three summands. A lexi¬ 
cographic order of the entries is again used, so that each row or column index k of 
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A=M®I®I + I®M®I + I®I®M corresponds to a triplet k = (A:i, ^ 2 , ^ 3 ) in 
the three-dimensional Cartesian grid. 

Corollary 6.3. Let M be (3-banded, Hermitian and with spectrum contained in 
[0,4p], and let A = M®I®I + I®M®I + I®I®M and k = {ki, ^2, fca) and 
t={tiA 2 -,H)- Then 

(exp(-T^))fet = (exp(-TM))fej,t,(exp(-rM))t2,/c2(exp(-TM))t3,fe3, 
from which it follows 


|(exp(-T^))fet| < $(fci,<i)$(fc 2 ,t 2 )$(fc 3 ,^ 3 ), 

for all {ki,ti), (* 2 ,^ 2 ), (ksAs) with min{|fci - <i|, \k 2 - < 2 !, 1^3 - ^ 3 !} > a/4tp/3. 
Proof. We write ^ = M®(/(g)/)-|-J(g) (M (g) J -|- J ® M), so that 

exp(—r^) = exp(— tM) 0 exp(—rM (Si I + I (Si {—tM)) 

= exp(— tM) C) exp(— tM) C) exp(— tM). 


Therefore, using z = (^ 1 ,^ 2 ), 


(exp(-T^))a = e^vec ((exp(-rM) 0 exp(-TM))eie^ exp(-TM)) 

= e^vec (vec(exp(—rM)et3e^ exp(—TM))e^ exp(—rM)) 


6^3 exp(-TM) 


= (efci exp(-TM)et3) exp{-TM)ek^) {ef^ exp(-TM)efe3). 



/ 

'exp{-TM)et,{el^ 

exp(—TM)ei) 

\ 

T 


exp(-TM)ei 3 (e^ 

exp(—TM)e 2 ) 


1 

_exp(-TM)ei 3 (ef 3 

exp(-TM)e„)_ 

/ 


The rest follows as in the proof of Theorem 16.II □ 

Remark 6.4. Using i 2 . S \} . one can obtain similar bounds for cos(^) and sin(^), 
where A = Mi ® I -\-1 ® M 2 with Mi, M 2 banded. 

6.2. Laplace—Stieltjes functions. If / is a Laplace-Stieltjes function, then 
f{A) is well-defined and exploiting the relation ( 12 . 21 ) we can write 


poo poo 

f{A)= / exp(—T.A)da(r) = / exp(—rM) 0 exp(—rM)da(r). 

Jo Jo 


Thus, using k = (fci,fc 2 ) and t = {ti,t 2 ), 


{f(A))kt = ei exp(-rM) ® exp(-rM)etda(T) 


0 

(exp(-TM))fc3ti(exp(-rM))t3fe3da(r). 

With the notation of Theorem 14.41 we have 



pOO 

|/(-4)|fct < / exp(-2AminT)| exp(-TM)|fe 3 t 3 | exp(-TM)|fe 3 t 3 da(T). (6.3) 
Jo 

In this form, the bound (16.3p . of course, is not particularly useful. Explicit bounds 
can be obtained, for specific examples of Laplace-Stieltjes functions, by evaluating or 
bounding the integral on the right-hand side of (j6.3l) . 
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For instance, using once again the inverse square root, so that a{T) = we 

obtain 

^ roo 2 ^ 

\A~^\kt<\/^ ^=exp(-2Amin'r)|exp(-TM)|fcitJexp(-rM)|fc2i2dr (6.4) 
Jo V'T'^ 

The two integrals can then be bounded as done in Theorem l4.4l For the other example 
we have considered earlier, namely the function f{x) = (1 — exp(—a;))/x, the bound is 
the same except that ^/n/\/ t^ is replaced by one, and the integration interval reduces 
to [0,1]; see also Example 16.51 next. 

Example 6.5. We consider again the function f{x) = (1 — exp(— a;))/a;, and 
the two choices of matrix M in Example 14.31 for each of them we build A as the 
Kronecker sum A = M®I + I® M. The entries of the tth column with t = 94, 
that is (^ 1 ,^ 2 ) = (14,5) are shown in EigurelH^l together with the bound obtained 
above. The oscillating pattern is well captured in both cases, with a particularly good 
accuracy also in terms of magnitude in the tridiagonal case. The lack of approximation 
near the diagonal reflects the condition |A:i — ti|/,5>2,i = l,2. 






Fig. 6.2. Examvle 1 6 . 51 True decay and estimates for \ f{A)\:^t, t = 94, A = M^I + I^M of 
size n = 400. Left: M = tridiagf— 1 , 4, — 1 ). Right: M = pentadiag(—0.5, — 1 , 4, — 1 , —0.5). 


Remark 6.6. These results can be easily extended to the case where A = Mi ® 
I + I ® M 2 with Ml, M 2 both Hermitian positive definite and having bandwidths fii, 
^2 ■ It can also be generalized to the case where A is the Kronecker sum of three or 
more banded matrices. 
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6.3. Cauchy—Stieltjes functions. If / is a Cauchy-Stieltjes function and A 
has no eigenvalues on the closed set F C C, then 


f{A) = j {A-ojT) ^ 7 ( 0 ;), 

so that 

elf{A)et = J^el{A-ujl)-^etdj{uj). 

We can write A—ujI = M®I+I®{M — u}I). Each column t of the matrix inverse, 
Xt := {col — A)~^et, may be viewed as the matrix solution Xt = Xt{uj) G to the 

following Sylvester matrix equation: 


MXt + Xt{M - wl) = Et, Xt=vec{Xt), et=vec{Et), 


where the only nonzero element of Et is in position (ti, ^ 2 ); here the same lexicographic 
order of the previous sections is used to identify t with (^ 1 ,^ 2 )- 

From now on, we assume that F = (—oo,0]. We observe that the Sylvester 
equation has a unique solution, since no eigenvalue of M can be an eigenvalue of 
00 1 — M for ct) < 0 (recall that M is Flermitian positive definite). 

Following Lancaster m p.556]), the solution matrix Xt can be written as 

pCO 

Xt = — exTp{—TM)Et exp(—r(M — a;/))dr. 

00 

For k = (fci, ^ 2 ) and t = (^ 1 ,^ 2 ) this gives 
el{ool - A)~^et = el^XtCk^ 

exp(-TM)etie^ exp(-T(M - w/))efe 2 dr. (6.5) 
Therefore, in terms of the original matrix function component, 

/ O noo 

/ e^i exp(-TM)et,ef^ exp(-T(M - ujI))ek 2 dTd'y{uj). 

-00 J 0 



We can thus bound each entry as 


|er/(.4)et| < 


exp(-TM)|fcitJ exp(-TM)|fc 2 t 2 / exp(Tw)d 7 (w) I dr 


It is thus apparent that \e'^f{A)et\ can be bounded in a way analogous to the case 
of Laplace-Stieltjes functions, once the term exp(Tw)d 7 (a;) is completely deter¬ 
mined. In particular, for f{x) = , we obtain 


fO 

/ exp(ra;)d 7 (a;) = 
J —00 
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Therefore, 


1 ^ \exp{-TM)\l^tJ{T)dT^ \exp{-TM)\l^tJ{T)dT^ (. 6 . 6 ) 

Using once again the bounds in Theorem 14.21 a final integral upper bound can be 
obtained, in the same spirit as for Laplace-Stieltjes functions. 

We explicitly mention that the solution matrix Xt could be alternatively written 
in terms of the resolvent (M — , with ^ G K |3S]. This would allow us to obtain 

an integral upper bound for \eJf{A)et\ by means of Proposition 14.71 and of (14.61) . 
We omit the quite technical computations, however the final results are qualitatively 
similar to those obtained above. 

Example 6.7. In Figure [6?3l we report the actual decay and our estimate fol¬ 
lowing (16.611 for the inverse square root, again using the two matrices of our previous 
examples. We observe that having used estimates for the exponential to handle the 
Kronecker form, the approximations are slightly less sharp than previously seen for 
Cauchy-Stieltjes functions. Nonetheless, the qualitative behavior is captured in both 
instances. 



Fig. 6.3. Example \ 6. 7| True decay and estimates for \A 2 |. f ^ t = 94 ^ A = M®I + I®M of 
size n = 400. Left: M = tridiag(—1, 4, —1). Right: M = pentadiag(—0.5, —1, 4, —1, —0.5). 


Remark 6.8. As before, the estimate for {f{A))k,t can be generalized to the sum 
A = Ml 0 / -|- 7 0 M 2 , with both Mi, M 2 Hermitian and positive definite matrices. 

Remark 6.9. Using the previous remark, the estimate for the matrix function 
entries can be generalized to matrices that are sums of several Kronecker products. 
For instance, if 


A = M®I®I + I®M®I + I®I®M, 


then we can write 

yl = M0(j0 7)-tJ0(M0/-tJ0M)=:M07-tJ0 M2, 
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SO that, following the same lines as in (16.51) we get 

elf{A)et = j^el{A- 

exp(-TM)“^etie^ exp(-T(M2 - ujI))ek2dTd'y{u}). 

Since M 2 = M ® I + I ® M, we then obtain exp(—TM 2 )efc 2 = exp(— tM) (g) 
exp(—TM)efe 2 . Splitting < 2,^2 in their one-dimensional indices, the available bounds 
can be employed to obtain a final integral estimate. 

7. Conclusions. In this paper we have obtained new decay bounds for the en¬ 
tries of certain analytic functions of banded and sparse matrices, and used these 
results to obtain bounds for functions of matrices that are Kronecker sums of banded 
(or sparse) matrices. The results apply to strictly completely monotonic functions and 
to Markov functions, which include a wide variety of functions arising in mathematical 
physics, numerical analysis, network science, and so forth. 

The new bounds are in many cases considerably sharper than previously published 
bounds and they are able to capture the oscillatory, non-monotonic decay behavior 
observed in the entries of f{A) when M is a Kronecker sum. Also, the bounds capture 
the superexponential decay behavior observed in the case of entire functions. 

A major difference with previous decay results is that the new bounds are given in 
integral form, therefore their use requires some work on the part of the user. If desired, 
these quantities can be further bounded for specific function choices. In practice, the 
integrals can be evaluated numerically to obtain explicit bounds on the quantities of 
interest. 

Although in this paper we have focused mostly on the Hermitian case, extensions 
to functions of more general matrices may be possible, as long as good estimates 
on the entries of the matrix exponential and resolvent are available. We leave the 
development of this idea for possible future work. 
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